function xyz=f_calibration(xy,pars)
CENTER_X=pars(1);
CENTER_Y=pars(2);
DX=pars(3);
DY=pars(4);
FOCALLENGTH=pars(5);
EP1=pars(6);
EP2=pars(7);
EQ1=pars(8);
EQ2=pars(9);
x=(xy(:,1)-CENTER_X)*DX;
y=(xy(:,2)-CENTER_Y)*DY;
r2=x.*x+y.*y;
dx=(x .* ((EQ1 * r2) + (EQ2 * r2 .* r2)))+((EP1 * (r2 + (2 * x .* x))) + (2 * EP2 * x .* y));
dy=(y .* ((EQ1 * r2) + (EQ2 * r2 .* r2)))+((EP2 * (r2 + (2 * y .* y))) + (2 * EP1 * x .* y));
xyz=[dx-x,dy-y,FOCALLENGTH*ones(size(x,1),1)];
xyz=xyz./repmat(sqrt(sum(xyz.*xyz,2)),1,3);

